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Although accumulation of molecular damage is suggested to be an important molecular mech¬ 
anism of aging, a quantitative link between the dynamics of damage accumulation and mortality 
of species has so far remained elusive. To address this question, we examine stability properties 
of a generic gene regulatory network (GRN) and demonstrate that many characteristics of aging 
and the associated population mortality rate emerge as inherent properties of the critical dynamics 
of gene regulation and metabolic levels. Based on the analysis of age-dependent changes in gene- 
expression and metabolic profiles in Drosophila melanogaster, we explicitly show that the underlying 
GRNs are nearly critical and inherently unstable. This instability manifests itself as aging in the 
form of distortion of gene expression and metabolic profiles with age, and causes the characteristic 
increase in mortality rate with age as described by a form of the Gompertz law. In addition, we 
explain late-life mortality deceleration observed at very late ages for large populations. We show 
that aging contains a stochastic component, related to accumulation of regulatory errors in tran¬ 
scription/translation/metabolic pathways due to imperfection of signaling cascades in the network 
and of responses to environmental factors. We also establish that there is a strong deterministic 
component, suggesting genetic control. Since mortality in humans, where it is characterized best, 
is strongly associated with the incidence of age-related diseases, our findings support the idea that 
aging is the driving force behind the development of chronic human diseases. 


Introduction 

Aging is a complex biological process. It occurs at 
every possible scale characterizing the living organism, 
ranging from the single cell level {e.g., oxidative damage, 
somatic mutations) to the level of interaction between dif¬ 
ferent organs {e.g., failure of individual organs with age 
or accumulation of dangerous byproducts of metabolic 
activity such as arterial plaque). This is why it has been 
notoriously difhcult to pinpoint the ultimate cause, a sin¬ 
gle molecular mechanism behind aging and relate it to 
such life-long consequences as the onset and progression 
of age-related diseases and, finally, death. As a conse¬ 
quence, a wide spectrum of propositions on the origin of 
aging has emerged over the years such as theories that 
view aging as a pre-programmed process [U l2 , various 
concepts involviM damage accumulation l3|-l7| , hyper- 
function theory |8l-[l0l| , disposable soma theory [TJ, , 
antagonistic pleiotropy |l3l | and mutation accumulation 
theory H. Since the theories of aging often differ greatly 
from each other both in spirit and letter and reflect upon 
different “facets” of aging iiEIil, it is not always 
possible to see their mutual inconsistencies and compat¬ 
ibilities. 

One of the key ideas in the study of aging is the causal 
role of molecular damage. Over 40 years ago Leslie Orgel 


proposed an attractive damage-based theory^, in which 
he argued that molecular corruption of enzymatic ma¬ 
chinery of transcription and/or translation would pro¬ 
duce positive feedback, leading to an “error avalanche”. 
We recently proposed a more general hypothesis of er¬ 
ror accumulationfT^ . supported by a semi-quantitative 
stability assessment for a generic gene regulatory net¬ 
work (GRN). We suggested a causal role of regulatory- 
error accumulation in the age-dependent increase in mor¬ 
tality rate, reflecting limited ability of cellular damage- 
response pathways to repair the consequences of inter¬ 
nal or external stresses However, no underlying mecha¬ 
nism has been established thus far, for the interplay be¬ 
tween the dynamics of error accumulation in transcrip¬ 
tion/translation/metabolic processes and the increase in 
mortality rates. 

To identify such a link, here we develop a theory of 
GRN stochastic critical dynamics and show that under 
very generic assumptions GRNs of most species should 
be inherently unstable in order to be compatible with 
age-dependent behaviour of mortality rates. Apply¬ 
ing the method of proper orthogonal decomposition [l8| 
to the publicly available age-dependent transcriptome 
and metabolome datasets of Drosophila melanogaster 
[T^ , we demonstrate explicitly that autocorrelation 

functions of transcriptional and metabolic datasets for 
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D. melanogaster exhibit stochastic exponential instabil¬ 
ity with a characteristic time scale, which is of the 
same order as mortality rate doubling time tMRDT of 
D. melanogaster. We provide explicit solutions for the 
stochastic aging dynamics of realistic GRNs to obtain ex¬ 
pressions for the age-dependent mortality increase. We 
argue that the two time scales, and tMRDT, should 
coincide for any species in which mortality follows the 
Gompertz equation. We are thus able to causally re¬ 
late GRN instability to the characteristic Gompertzian 
increase of the mortality rate in populations. The expo¬ 
nential increase of mortality eventually slows down and 
the mortality rate is expected to approach a plateau, 
^ ^is) 1/iMRDT, at late ages. The result should 
be valid for sufficiently long lived species, t\s » ta- 
The prediction is conhrmed by quantative analysis of 
the transcriptomes, metabolomes and mortality curves 
for populations of D. melanogaster and mortality curves 
for very large cohorts of medflies [2l|. We establish spe¬ 
cific metabolites and genes as biomarkers of aging in D. 
melanogaster., examine the respective human orthologs, 
and provide their associations with genes and pathways 
commonly related to the major human chronic diseases. 

We show that the dynamics of aging, as mani¬ 
fested in transcriptional and metabolic profiles, can 
be decomposed into both stochastic component, re¬ 
lated to regulatory error accumulation in transcrip¬ 
tion/translation/metabolic pathways, and strong deter¬ 
ministic component, which can be naturally associated 
with a genetic program of an organism. We show that 
the presented model may serve to embrace most of the 
known features of the aging processes and hence provides 
a novel and universal basis for quantitative analysis of di¬ 
verse age-dependent processes. We believe that the theo¬ 
retical ideas presented in the paper could be useful in the 
identification of biomarkers and, after appropriate devel¬ 
opment, in mechanistic studies of processes that regulate 
aging in any sufficiently long-lived organism. 

RESULTS 

Critical dynamics of Gene Regulatory Networks: 

quantitative model of aging 

In this Section, we first explain the theoretical basis of 
our quantitative analysis of transcriptional and metabolic 
profile changes with age and construct a quantitative 
model of aging. We consider the case of the GRN de¬ 
scribed by a generic set of non-linear matrix equations of 
systems biology 

g [x,dx/dt,(fx/dt^,...) = F. 

Here 5 is a vector function of any number of quanti¬ 
ties representing the state of the system, x, and its time 
derivatives which characterize the interactions between 
different components of x. If, for example, the state vec¬ 
tor X consists of the gene expression levels only, e.g. as in 


, the function g may be thought of as encoding path¬ 
ways to which the given genes contribute. In addition 
to gene expression levels, the state vector x may include 
other descriptors, such as levels of metabolites or 
methylation levels in DNA [2^ . As the system state vec¬ 
tor X is not observable in its entirety, we presume that 
any sufficiently large part of it captures representative 
information about “macroscopic” properties of the bio¬ 
logical system state at the organismal level. For the sake 
of simplicity, below we will refer to gene expression only, 
if not stated otherwise. 

The vector F describes the action of external or inter¬ 
nal stress factors affecting components of x. The vector 
F may depend on time, t: F(t) = Fq + where Fq 

and SF represent the mean stress and the fluctuations of 
stress levels, respectively. 

Over long times, gene expression levels typically fluc¬ 
tuate near certain, relatively slowly changing, mean val¬ 
ues, corresponding to the average homeostatic state of 
an organism. Therefore, we assume that there is always 
a quasi-stationary point xq (homeostasis of the organ¬ 
ism), given by the solution of the stationary equation 
/(cco) = Fq- As a result, the slow dynamics of the fluc¬ 
tuations of the gene expression levels, Sx = x — Xq, in its 
leading order obeys the matrix stochastic equation for a 
linear continuum-limit of a Markov chain: 

DSx + KSx = SF, (1) 

where the matrices D and K describe the dynamical re¬ 
laxation properties and the interactions between the com¬ 
ponents of the gene regulatory network, respectively (see 
Supplementary Information, Section for in-depth dis¬ 
cussion of Eq. (ID) derivation). 

It has been previously suggested that, in most species, 
GRNs operate close to a stability-instability transition 
or order-disorder bifurcation (TtI . I^ . . The transition 

from stability to instability in networks with the network 
graphs not possessing specific symmetries is typically as¬ 
sociated with occurrence of co-dimension I bifurcations 
[2^ [2^ . Such transitions are characterized by the loss 
of stability along a single direction in the state vector 
space, corresponding to the first principal component, 
while contributions from all other principal components 
remain stable. This situation, known in the literature 
as a saddle-node bifurcation [^, is realized when the 
smallest real eigenvalue, e, of the matrix K tends to zero, 
e —?► 0, and then becomes negative. Accordingly, the sta¬ 
tionary solution for xq ceases to exist, and the homeo¬ 
static state of the organism starts to change in time. As 
we shall argue, this very instability is directly responsible 
for the process of aging. 

To be specific, hereinafter we focus on this scenario. 
As explained in Supplementary Information, Appendix 
El the derivative 5x of any gene expression level becomes 
critical close to the co-dimension 1 bifurcation point. 
This implies that fluctuations of the expressome 5x are 
strongly amplihed with age, the phenomenon known as 
critical slowing-down (see e.g. |27l|'). In particular, the 
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Figure 1: (A) Schematic representation of the stochastic dynamics in the effective potential V{z) corresponding to 
the cases of unstable {a > 0) and stable (a < 0) GRNs, see Eq.©. (B) Age dependence of zi 2 = • 61 ^ 2 ) of the 

expressome state projections onto the vectors 61^2 corresponding to the first two principle components as functions 
of age (gene expressions for normally fed D. melanogaster from [H, red, and metabolite levels for normally fed D. 
melanogaster from blue). Here E{...) stands for the average over the biological repeats; (a^ • b) denotes a scalar 

product of two vectors a and b. 


autocorrelation time scale Guto = ■ Db)/e of the ex¬ 

pressome diverges as e —>■ 0 , indicating that stochastic 
dynamics of 5x is very slow near the point of bifurca¬ 
tion. The fluctuations of the expressome state vector 5x 
including those created in response to persistent external 
stress, are mostly amplified along the direction b of the 
right eigenvector of matrix K corresponding to the van¬ 
ishing eigenvalue e. The latter observation is supported 
by the analysis of experiments in D. melanogaster^ where 
transcriptional responses to different stress factors were 
found to contain a good fraction of similar differentially 
expressed genes in the limit of both nearly lethal [ 2 ^ and 
weak stressors (^ . 

Accordingly, over times comparable to the average 
lifespan, the gene expression fluctuations are dominated 
by transcript-level changes along the vector correspond¬ 
ing to the leading principal component, which is collinear 
with the direction of the vector b. We conclude that for 
long time intervals evolution of the vector of the GRN 
state vector, dx, can be accurately described using a sin¬ 
gle variable z, such that 5x k. z ■ b. Accordingly, Eq.([T]) 
can be reduced to the equation 


dz dV{z) 
dt dz 


( 2 ) 


describing the Brownian motion of a particle in an ef¬ 
fective potential V{z) ~ —a 2 :^/ 2 , for sufficiently small 
z, and the quantity a = —ejlaF ■ Db) characterizes the 
“stiffness” of the gene regulatory network. The stochastic 
force, / = (a^ • SF)/{a^ ■ Db), is a single unihed quan¬ 
titative measure of stochastic component of all external 
and internal stress, and is characterized by the diffusion 
coefficient A (see Supplementary Information, Appendix 
l^for the in-depth derivation and necessary justification). 

Depending on the specific morphological properties of 
the gene regulatory network, the curvature of V{z) at 
small 2 ; can be negative, a < 0 , or positive, a > 0 , as 
shown in Figure [T]4. The first case will be briefly consid¬ 
ered in Discussion section. In the latter case, the GRN 
is inherently unstable, and the variance of gene expres¬ 
sion levels, C{t), computed using the solutions of Eq.Q, 
grows exponentially with age as 


C{t) = E{Sx{t)Sx^{t)) ~ E{SxQ)b ■ 6 ^ exp(2at). (3) 

The expression ([3]) is fully consistent with earlier observa¬ 
tion of increasing biological variability with age (13, [3l|. 

Griticality of gene expression-level dynamics criticality 
implies that the parameter a is close to zero and hence 
the expressome covariance matrix C{t) should be a ma- 
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Figure 2: Basic biological evidence for the presented stochastic model of aging: (A) Mortality rate along with the 
projection z of the system state vector, Sx, represented by the transcriptomes of aging D. melanogaster 0 , onto 
the singular direction b of the covariance matrix E{5x{t)5x'^{t')) as a function of age. (B) Components of the vector 
b for D. melanogaster. Peaks correspond to the values of the components of the vector b for every gene in the 
dataset. (C) Mortality rate and the projection z of the system state vector 6x, represented by the metabolomes of 
D. melanogaster as a function of age (20|; (D) Components of the vector b for D. melanogaster, computed using only 

the levels of targeted metabolites 0. 
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Figure 3: Biological characterization of the “leading direction of aging” b in D. melanogaster: (A) Gene Ontologies 
(Biological Process only) enrichment by the leading positive and negative components of the vector b computed 
using gene expression transcripts; (B) KEGG pathways enrichment for human orthologs of the gene sets 
representing the leading components of b, corresponding to positively and negatively regulated genes. 


trix with special structure. Eq.Q implies that in this 
case the right eigenvector of the covariance matrix corre¬ 
sponding to its largest eigenvalue coincide with the right 
eigenvector b of the matrix K corresponding to vanishing 
eigenvalue e —> 0 (see Supplementary Information, Ap- 
pendix[A|. What is more, the largest eigenvalue should be 
much larger than all other eigenvalues of this matrix. In 
reality the properties of the experimental sample covari¬ 
ance matrix may be degraded by the experimental noise 
and necessarily very limited number of samples compared 
to the number of transcripts or metabolites observed. 

These conclusions are supported by Principal Gom- 
ponent Analysis (PGA) of the age-dependent shifts of 
RNA-transcript and metabolite levels in normally fed D. 
melanogaster from (T^ and , respectively. The covari¬ 
ance matrices computed from such datasets are indeed 
almost singular, see Supplementary Information, Figure 
m Correspondingly, the variation of the state vector Sx 
with age is dominated by the change along the loading 
vector corresponding to the first principal component, 
and z = zi, where zi = {Sx'^ ■ b) is the hrst princi¬ 
pal component, as described in Figure [Tjl. The variance 
along this vector grows very quickly with age, both in the 
metabolome and the transcriptome, in accordance with 
Eq.(l3]). The exponent, roughly inferred from each of the 
datasets for ages less than the mean lifespan, appears to 
be nearly the same, which is an indication of the same 
GRN instability observed in each of the experiments. For 


the sake of comparison, we have also plotted projections 
of gene expression on the next principal component. It is 
easy to see that the transcriptome and metabolome re¬ 
main dynamically stable along the direction of the second 
principal component. 

We thus hypothesize that at the origin, z ^ 0, the 
expressome state corresponds to a healthy or “youthful” 
state, while larger values of z phenotypically describe dis¬ 
torted states of aged animals. Accordingly, aging mani¬ 
fests itself as a slow dynamics of expression levels, an ex¬ 
ponential roll-out along the singular, or “aging” direction 
6 , from the younger animals to the older ones. The sin¬ 
gular “direction of aging” b is associated with a response 
to a generic stress and determined by the properties of 
the interactions in the network only. At least within the 
model assumptions, the direction b does not depend on 
the detrimental stress factors and appears to be hard¬ 
coded in the genome, and hence the development along b 
can be considered as an organism-level manifestation of 
an aging program or quasi-program The “direc¬ 

tion of aging”, vector 6, is similar but not the same as the 
differential expression vector in aging animals, since the 
latter may include a contribution of faster modes repre¬ 
sented by other eigenvectors of the GRN connectivity ma¬ 
trix K corresponding to eigenvalues with larger absolute 
value. Since the network state fluctuations along all the 
other directions are not related to aging and are relatively 
small, we conclude that the deterministic component of 
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aging in Sx(t) already dominates over the stochastic one, 
even early in life (also see below). 

Gompertz mortality law and biological age 

To relate the GRN instability to the process of aging 
we further assume that for every gene in the genome there 
exists a particular threshold associated with its over- or 
under-expression, which leads to a strong build-up of 
deleterious changes, which in it turn ultimately lead to 
death. Since fluctuation of the expression levels can be 
written in the form 5x k, z ■ b, such deleterious thresh¬ 
olds can be defined by the specific values oi z = Zi for 
every gene, which may be either positive or negative. It 
is reasonable to assume that the death of a cell (and ulti¬ 
mately of an organism) occurs once z reaches the small¬ 
est (by absolute value) of the threshold values, z = Z. 
In Supplementary Information Section |B] we argue, that 
as the expressome state deviates from the initial point, 
the higher order non-linearities in Eq. © and therefore 
in the effective potential V(z) of Eq.® can no longer be 
neglected. As expected from the argument, we are able 
to show in Supplementary Information section [b 1 that 
the threshold value Z can be equally well introduced for¬ 
mally in the model as the single effective property of the 
non-linearity in such a way that Z is large whenever the 
non-linear interactions are weak.. 

Initially, the probability density along the aging direc¬ 
tion is localized at a small value of z and the mortality is 
also very small. If severe deleterious changes occur late in 
life, i.e. if 7 = aZ^/A ^ 1 , the lifespan greatly exceeds 
ta, the non-linearity is small, and the mortality increases 
exponentially with age first, then the growth decelerates, 
and, eventually, the mortality reaches a plateau. At in¬ 
termediate ages, in a narrow range close to the average 
lifespan of the animals, the mortality can be approxi¬ 
mated with 

M « (4) 

This is a form of the Gompertz law (see Supplemen¬ 
tary Information, Appendix |B] for the complete analyti¬ 
cal solution for the age-dependent mortality in a closed 
form for all ages). The Gompertz exponent a is related 
to the Mortality Rate Doubling Time (MRDT), tMRDT, 
by equation a. = log(2)/tMRDT- In practice, the value 
of a is often determined from the survival records at the 
age, corresponding to a narrow interval of width ^ 1 /a 
around the mean lifespan, when most of the animals die. 
Therefore Eq.(|3]) suggests that the empirical Gompertz 
exponent should be very close to the GRN “stiffness”. 

On the other hand, the Initial Mortality Rate (IMR), 
Aimr, is related both to the stress levels and to the time 
required to reach the boundary Z under the action of the 
random damaging forces only, Aimr = ^jZ"^. Normally, 
the stochastic forces are weak, 7 ~ (o/Aimr) is large, 
and, therefore, average lifespan, E{t\s) ~ a“^log 7 ^ 


a“^, is also large compared to the characteristic time, 
associated with the gene network instability, ta ~ a~^. 
Remarkably, the lifespan in the model depends on the 
morphological properties of the GRN, such as its topol¬ 
ogy and connectivity, through the matrix AT, and, more 
specihcally, through its lowest eigenvalue e ^ a. The de¬ 
pendence of species lifespan on damaging factors through 
the parameter A is, on the contrary, logarithmically 
weak. This has been recently shown in analysis of data 
from [33 in our work[l3, where long-term effects of trau¬ 
matic brain injury at young age led to significant, but 
small changes in organismal lifespan. This also means 
that the gene-network evolution along the aging direc¬ 
tion is practically deterministic in sufficiently long-lived 
organisms. At the same time, the stochastic forces are 
still very important early in life, yielding a relatively wide 
distribution of transcription profiles even in a group of 
otherwise identical animals at birth. Our analysis pro¬ 
poses that the stochastic component of the GRN dynam¬ 
ics explains the distribution of the observed ages at death 
in such a cohort. 

Mortality in a form similar to Eq. (|3]) can be associated 
with gene network instability in a simple phenomenologi¬ 
cal model, where both IMR and MRDT were expressed in 
terms of generic network parameters, such as translation, 
gene repair and protein turnover rates, along with the 
stress levels, whereas mortality is directly associated with 
the increasing number of regulatory errors [13. There¬ 
fore, the quantitative indicator of the expressome remod¬ 
eling along the aging direction, z, can be interpreted as 
a measure of damage accumulated over the lifespan of an 
organism. 

As we show in Supplementary Information, Appendix 
IB the dynamics of gene expression levels 5x in the weak 
non-linearity, or “Gompertzian” ,limit, 7 ^ 1 , proceeds 
along a well-defined trajectory, corresponding to a con¬ 
tinuation of the development program. The “distance 
traveled” along the aging direction up to current mo¬ 
ment of time, z ^ {5x'^ ■ b), is, according to Eq.(I31), di¬ 
rectly related to mortality and hence by itself is a good 
biomarker of mortality and aging. These observations 
pave the way to define biological “clocks” using any num¬ 
ber of variables characterizing a GRN state. For exam¬ 
ple, Figures [5J4 & C suggest the possibility of developing 
transcriptome- and metabolome-derived biological clocks 
for D. melanogaster. Those clocks in turn may be com¬ 
pared with the biomarkers of aging and the biological 
age calculator (biological clock) similar to, for example, 
the one derived from a regression model established using 
age-dependent DNA-methylation patterns 


Association between GRN instability and mortality 

rate 

Our analysis of D. melanogaster age-dependent tran¬ 
scriptional data from [l^ indicates that the gene ex¬ 
pression variance grows exponentially with age in agree- 
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ment with Eq.(|31), with a k. 19.1 yr“^, as shown in Fig¬ 
ure [214. In a similar fashion we have studied the “direc¬ 
tion of aging” b in the metabolome of D. melanogaster 
[ 23 . The projection of the metabolome state vector onto 
the singular direction b increases very quickly with age, 
along with the reported mortality as shown in Figure 
WP. The mean lifespan of D. melanogaster is about 25 
and 40 days in both experiments, respectively, and in 
the two cases the initial exponential growth of the pro¬ 
jection z‘^{t) ceases at t > tis, as predicted by the the¬ 
ory. Figure [23 represents a side-by-side comparison of 
the variance computed from the gene expression and the 
metabolome datasets. The values of instability rates, a, 
recovered from analysis of both experiments, are very 
similar and close to the value of the Gompertz exponent 
a = ln2/tMRDT = 17.1 yr“^ corresponding to the mor¬ 
tality rate doubling time tiuRDT ~ 0.04 yr [Sg. 


The pace of mortality-rate increase slows at 
advanced ages 

At more advanced ages, t > E{t\s), where t\s is the 
average lifespan of the population, solutions of Eq.([2]) 
exhibit deceleration of mortality. More specifically, as 
shown in Supplementary Information, Appendix m we 
expect that the mortality rate ceases growing and satu¬ 
rates at a constant value, M^o ot. This quantitative 
prediction can be compared with the experimental data 
for medflies [2lj, where the cohort sizes were sufficient 
to observe both the regime of exponentially increasing 
mortality and the mortality plateau (see Supplementary 
Information, Figure Ej). For the male medflies we find 
the asymptotic value Moo ~ 0.16 day“^, which is close to 
the Gompertz exponent a ~ 0.22day~^, in accordance 
with theoretical prediction. The asymptotic value of the 
mortality rate for female medflies. Moo ~ 0.12 day”is 
also of the right order of magnitude, though nearly two 
times smaller than the slope of the mortality exponent 
a = 0.25day”^. Nevertheless, there is an approximate 
agreement between between the asymptotic value Moa of 
the mortality rate at late ages and the value of the Gom¬ 
pertz exponent a in both cases. We also note that the 
approximate Eq.([3|) should work better for longer living, 
or “Gompertzian”, Aimr a, organisms. 


Biomarkers of aging in D. melanogaster 

We computed the aging direction based on the tran- 
scriptomic data from Drosophila melanogaster 0- Gene 
Ontolology (GO) analysis of the leading components of 
b vector is summarized in Figure [314, for Biological Pro¬ 
cess (BP) categories only. The genes providing the most 
significant negative contribution to the aging direction b 
can be roughly split into three groups. First, there is 
a large group of genes related to metabolic processes: 
organophosphate metabolism, generation of precursor 


metabolites and energy, membrane and mitochondrial 
ion transport. A detailed review of the GO categories in¬ 
volved seems to indicate that most of the down-regulated 
processes are related to oxidative phosphorylation or res¬ 
piration. Another large group of genes involves devel¬ 
opmental processes, such as anatomical structure forma¬ 
tion involved in morphogenesis and extracellular struc¬ 
ture organization. Genes encoding eggshell and egg coat 
formation correspond to strongly negative components of 
the vector b, which points to the reproductive senescence 
in female flies. Finally, we observe decline in processes 
related to proteostasis, including protein folding, ER lo¬ 
calization and proteolysis. 

As expected, the cluster of genes corresponding to the 
strongest positive components of the vector b includes 
genes associated with a generic stress response. Among 
other leading components there are genes corresponding 
to the innate immune response (which is also a stress 
response), in agreement with previous conclusions ( 3 ^ . 
and the components of oxidation-reduction and reactive 
oxygen species metabolic processes. 

We found that several transcription factors (TF) ap¬ 
pear to be strongly associated with aging. One of the 
top positively associated genes, hairy, is a transcrip¬ 
tional suppressor, and believed to be a metabolic switch, 
found up-regulated in the gene expression datasets char¬ 
acterizing oxygen-deprived flies. Gonversely, mutations 
in hairy significantly reduce hypoxia tolerance (^ . The 
only negatively associated transcription factor, eve, is a 
transcription factor involved in multiple organ and sys¬ 
tem development. 


Human orthologs for biomarkers of aging in D. 
melanogaster 

To obtain a better insight into aging of D. melanogaster 
we annotated gene expression changes involved in aging 
in flies with human orthologs. This is possible to do 
since according to OrthoDB, approximately 80% of all 
genes in D. melanogaster have human orthologs. We 
searched for human orthologs of the leading contribu¬ 
tors to the aging direction vector b computed from the 
transcriptome of D. melanogaster. The results of this 
analysis are represented in Figure[3l3. Remarkably we ob¬ 
serve transcriptional changes associated with human age- 
related diseases: the downregulated genes were enriched 
in KEGG categories corresponding to age-related neu- 
rodegenerative diseases (Huntington’s, Parkinson’s and 
Alzheimer’s), and also the processes of cardiac muscle 
contraction. We also observed that the downregulated 
genes enriched in KEGG pathways were associated with 
oxidative phosphorylation (consistent with the GO an¬ 
notation and KEGG enrichment results above). In con¬ 
trast, KEGG PPAR signaling pathway is activated with 
age, which may be another indication of the activation of 
anti-oxidant metabolism systems as noted above. There 
is evidence suggesting that PPARa and the genes under 


its control play a role in the evolution of oxidative stress 
increases observed in the course of aging in mice (dlj . 

Further evidence relating the stochastic instability of 
the gene regulatory network to the common age-related 
diseases is a strong positive contribution of the Pepck 
gene, which encodes a rate-controlling enzyme of glu- 
coneogenesis, the process by which cells synthesize glu¬ 
cose from metabolic precursors. RNAi silencing of the 
Pepck gene was found to be an effective way to reverse 
diabetes-induced hyperglycemia in miceji^. It has also 
been shown that expression of the Pepck ortholog in 
Caenorhabditis elegans correlates almost perfectly with 
longevity in an isogenic series of longevity mutants . 
Also, a transgenic mouse strain constitutively overex¬ 
pressing muscle PEPCK is muscular and long-lived (43 |. 


DISCUSSION 

Our analysis suggests that aging involves an organis- 
nial level manifestation of the inherent instability of gene 
regulatory networks. Proper orthogonal decomposition 
of age-dependent transcriptional and metabolite prohles 
shows that, for species which age in accordance with the 
Gompertz equation, the life-long dynamics of the tran¬ 
scriptional and metabolite profiles can be effectively de¬ 
scribed in terms of stochastic critical dynamics with a 
single degree of freedom; that dimension/component is 
identified by the projection of the expressome on the 
right eigenvector b of the CRN connectivity matrix K, 
and associated with the eigenvalue e which vanishes at 
the critical point. We establish that depending on the 
sign of e, a CRN can be either stable or unstable. We fur¬ 
ther propose that changes in metabolite levels contribute 
to aging similarly to changes in CRN, e.g., through a 
buildup of certain metabolites and molecular damage be¬ 
yond a certain level. We note that appearance of such 
deleterious thresholds is guaranteed, if non-linearities in 
the equation o for expressome levels are accounted for 
(see Supplementary Information, Appendix |B|) . We show 
that in this limit mortality follows the form of the Gom¬ 
pertz law, i.e. it first increases exponentially with age 
and then saturates at a constant level, a prediction sup¬ 
ported by experimental evidence. As predicted, we hnd 
approximate agreement between the rate of gene regula¬ 
tory network instability and the value of mortality rate 
saturation. Thus we are able to demonstrate that Gom- 
pertzian aging is a property of organisms with inherently 
unstable gene regulatory networks, and thus the gene net¬ 
work instability is a process that directly contributes to 
aging, which can be observed in the form of an exponen¬ 
tial increase of all-cause mortality with age. 

Though stochasticity of the expressome is an essen¬ 
tial part of the model, it should be noted that in the 
Gompertzian limit a wide variety of organisms exhibits 
similar dynamics of expression prohles. According to our 
theory, genetically programmed connectivity matrix K 
can be represented by its low-rank approximation using 


its right eigenvector b and the corresponding eigenvalue e. 
Such a direction can be interpreted as a genetically pro¬ 
grammed mechanism of aging. It is important to note 
that this direction is identihable from experimental data, 
thus allowing quantitative predictions to be made about 
a mechanism behind aging. 

The stochastic model of aging presented here is fully 
compatible with, and in fact embraces, several prevailing 
theories and hypotheses of aging, in the following sense: 

• Prc^rammed and quasi-programmed aging 
EiO- We have shown that aging-driven dy¬ 
namics of the expressome 6x(t) can be separated 
into stochastic and deterministic components. The 
deterministic component of Sx(t) starts to domi¬ 
nate over the stochastic one even early in life, at 
roughly the same time as the low-rank approxima¬ 
tion of the GRN connectivity matrix K becomes 
prominent, as can be seen from the fact that the 
covariance matrix is singular and dominated by the 
aging direction. As the dynamics of the expressome 
Sx(t) is largely determined by the properties of the 
matrix K, we interpret the low rank approximation 
of K as a program or a quasi-program of aging Q , 
though no statement about its evolutionary nature 
can be made. We speculate that the same direction 
that corresponds to aging occurs and is active dur¬ 
ing development, thus being consistent with quasi- 
programmed aging theory. 

• hyperfunction theory [^. The single mode 
6x = z ■ b, describing a large cluster of genes that 
make a dominant contribution to the process of ag¬ 
ing, starts dominating the aging dynamics soon af¬ 
ter the completion of development, and the stochas¬ 
tic effects on the expressome dynamics can be con¬ 
sidered relatively weak i,a. The expression of 
positive leaders of the aging direction b will grow 
with time, which can be interpreted as hyperfunc- 
tion. In the model proposed, death of the organism 
is associated with a threshold beyond which this or¬ 
ganism dies, which can be interpreted in the man¬ 
ner of hyperfunction theory as excessive functions 
of gene products for positive leaders of aging di¬ 
rection. However, some signihcant differences from 
hyperfunction theory as proposed in ii should 
be noted. Namely, we do not have any experimen¬ 
tal data to determine the role of aging direction in 
developmental processes so this connection remains 
speculative. Moreover, according to the theoretical 
model, changes along stable developmental direc¬ 
tions will tend to diminish with time and will not 
contribute to changes in gene expression, so at least 
some fraction of genes involved in development will 
not contribute to aging. 

• damage/error accumulation (^-Q. The 

stochastic component of the expressome Sx(t) can 
be interpreted as a result of accumulation of er¬ 
rors, e.g. regulatory errors. Though it becomes 
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small compared to the deterministic component at 
late ages, its effects are strongly pronounced at all 
ages, as the mean square deviation yjE[5x^{t)) first 
grows with age at the same exponential Gompertz 
rate as the mortality M{t) and later stabilizes at 
some constant level, contributing at late ages to the 
very process of organism passing beyond toxicity 
threshold and dying. Our analysis shows, that the 
stochastic forces shaping the expressome evolution 
early in life explain a wide distribution of the ages 
at death even in a genetically homogenous popula¬ 
tion. It should be noted that damage accumulation 
does not exclude the scenario of quasi-programmed 
aging, as damage is not limited to byproducts, er¬ 
rors and other molecular forms, but encompasses 
deleterious changes at all levels, or the deleteriome 


The relationship just established between the critical dy¬ 
namics of the gene regulatory network state vector and 
the Gompertzian mortality characteristic of most species 
allows one to consider the projection z = • 5x) as a 

candidate biomarker of aging. Due to the low value of the 
eigenvalue e the direction h dominates the response to any 
generic stress. This, in fact, was already demonstrated 
in (^ . where transcriptional signatures of responses to 
very different stresses were found to share a great number 
of gene expression changes. 

The levels of gene expression or metabolites corre¬ 
sponding to the most significant components of the ag¬ 
ing direction b are strongly associated with age and age- 
related diseases. Therefore aging as described here, in 
terms of GRN instability, leads to the impairment of nor¬ 
mal physiological functions, has characteristic biomark¬ 
ers and phenotypes including signs of major diseases of 
aging, and produces exponentially increasing morbid¬ 
ity. Accordingly, the process of aging itself falls un¬ 
der the definition of a disease recently used by AMA 
for a common condition such as obesity [i^. We must 
note here that therapeutic or experimental interventions 
aimed to counter-balance the age-related changes in gene 
or metabolite expression may not be very effective ways 
to extend lifespan of the species. For example, even 
though aging in D. melanogaster is associated with in¬ 
creased internal and external bacterial load, as well as 
with increased expression of antibacterial peptides, nei¬ 
ther reducing the bacterial population with antibiotics 
nor reducing the humoral antibacterial response was able 


to extend lifespan j^. The transcription factor hairy is 
overexpressed with age, a nd y et its inhibition does not 
result in lifespan increase [^. This reinforces the con¬ 
clusion that the markers of aging are, in general, not the 
same as regulators of aging. 

Even though most of the analysis in this manuscript 
is performed for D. melanogaster^ our conclusions are 
generic and should be applicable to other species. De¬ 
pending on the network and environmental parameters, 
realistic GRNs would not necessarily support a deriva¬ 
tion of the Gompertz law. In fact, aging is extremely 
diverse across the tree of life [d^. We believe that one of 
the most intriguing form of a non-Gompertzian mortality 
law in the model may arise if the effective potential V(z) 
has a local minimum with small but positive curvature, 
as in Figure [TJ4, a < 0. The higher order nonlinearities, 
such as the cubic terms in the effective potential, cannot 
be neglected anymore and in this case the genetic network 
turns out to be metastable. If the minimum of the poten¬ 
tial R(z) is separated from the region of large z by a suffi¬ 
ciently high activation barrier (see Figure[T]4 and Supple¬ 
mentary Information, Appendix |B]) , the mortality rate, 
determined by the probability of activation, is exponen¬ 
tially small and age-independent. This situation could be 
a physical picture behind the hypothesis that some ani¬ 
mal species exhibit negligible senescence [d^ This feature 
of our model may also explain the apparent lack of age- 
dependent changes in physiological parameters as well as 
mortality rates observed in the naked mole rat and some 
other species over long time periods [d^. This argument 
may also be supported by the analysis in nam , where 
the number of the genes differentially expressed with age 
was compared between extremely long-lived species and 
mice and humans. These works show that the number 
of differentially expressed genes in species with negligible 
senescence is much lower than in mammals with Gom¬ 
pertzian aging. 
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SUPPLEMENTARY INFORMATION 

Appendix A: Stochastic dynamics of gene expression 
and metabolite levels 

In the following, we explain the physical basis for Eq. 
© and provide its derivation. We also explicitly find 
the expressions for the mortality rate, the survival prob¬ 
ability, the lifespan distribution function and the average 
lifespan of species with expressome subject to Eq. 

1. Getting a grasp on expressome dynamics. 
Model reduction and physical considerations 



The analysis represented in this paper is mainly based 
on the observation that time series datasets describ¬ 
ing changes of expressome profiles with age are of¬ 
ten quite susceptible to Model Reduction techniques 
[T^ . Proper orthogonal decomposition of dynamic (time- 
series) datasets for age-associated gene expression lev¬ 
els (and also metabolite levels or other “omics” datasets) 
[l^ 1^. [ 5 ^ shows that the long-time behavior of any ex¬ 
pressome X is largely determined by a single component 
(see Fig. H). The number of different principal compo¬ 
nents PCn representing the signal x(ti),... can¬ 

not not exceed the number of time snapshot points in 
corresponding time series. Even though m is typically 
limited and is very small compared to the number of 
transcripts or metabolites, , the clear dominance of a sin¬ 
gle component in the signal is a non-trivial fact. As we 
argue here, it carries important information about bio¬ 
logical gene regulatory and metabolic networks; namely, 
it shows that such networks are pre-critical in a sense 
that will be explained below. 

Consider a meta-stable gene regulatory network 
(CRN) with the instantaneous network state represented 
by a vector x, whose components Xi are given by expres¬ 
sion levels of different genes, proteins or metabolites par¬ 
ticipating in the CRN. Dynamics of the vector x consists 
of mRNA levels, concentrations of proteins and metabo¬ 
lites, involved in regulatory pathways and are influenced 
by external as well as internal stress factors acting on 
the cell. This dynamics is governed by differential ma¬ 
trix equations of systems biology 


/ dx d^x \ 


= F, 


(Al) 


where g is a vector function of the state vector x, en¬ 
coding interactions between different components of the 
expressome, and the vector F describes the action of 
(mostly stochastic) external or internal stress factors af¬ 
fecting components of x. For the sake of simplicity, below 
we will refer to gene expressions only, if not stated oth¬ 
erwise. 

Generally speaking, the vector function g may depend 
on arbitrarily high time derivatives ^^,..., ,... 


Figure 4: Five largest singular values Si{x{t)) of the 
expressome x(t) divided by the largest one, Si(t)] 
transcriptome of D. melanogaster 0 (blue), 
transcriptome of R. norvegicus, hippocampus |52l | 
(green), metabolome of D. melanogaster, targeted 
metabolites only 0 (red), metabolome of D. 
melanogaster, all metabolites [20| (purple). One of the 
orthogonal components of the expressome x{t) clearly 
dominates over the others in all considered cases. 


of the expressome state vector. Hereinafter the focus of 
the study is on the dynamics of x at time scales compara¬ 
ble to the mortality rate doubling time tiviRDT; which is 
of the order of tens of days for nematodes and drosophi- 
lae and hundreds of days for mice. This time scale is to 
be compared to a typical time scale of protein transla¬ 
tion and regulation, or a lifetime of mRNA. In order to 
describe this slow dynamics of the expressome state vec¬ 
tor X, one can neglect all time derivatives of x higher the 
first derivative ^. The matrix Eq. thus reduces to 

(A2, 


2. Prequency domain properties of the vector F of 
stress factors 

As it was explained in the main text, the vector F of 
environmental and internal stress factors affecting gene 
expression levels x has both a slowly changing compo¬ 
nent, Fb(0) ^nd a stochastic component, SF(t), rapidly 
changing in time. Below we assume for simplicity that 
Fo(0 = Fq, i.e., the slowly changing component is, in 
fact, constant. 

Since dynamics of the vector SF can be considered es¬ 
sentially random at time scales of interest, dynamics of 
X is driven by correlation properties of SF. The total 
number of environmental and internal stress factors af¬ 
fecting gene expression levels x is extremely large, and 
the vector SF{t) represents a superposition of these fac¬ 
tors. According to the central limit theorem it can 
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thus be considered a Gaussian stochastic process, so that 
E{6F{t)SF{t')) = B{t-t'), 

where B{t — t') is a function of the time difference \t —1'\ 
only, while F{...) stands for the expected value of its ar¬ 
gument. The statistical average is understood as average 
over ensemble of species/cells of a given organism. In the 
frequency domain, one has correspondingly 

E{5FU)5F{-f)) = 

= J df e-‘^^^^<'*-*'^F{6F{t)SF{t')) = B{f), 

where / is frequency of a given mode in the Fourier ex¬ 
pansion of the function B{t — t'). 

The function B{f) has multiple singularities in the 
complex plane of /. These singularities encode charac¬ 
teristic time scales of pathways’ dynamics. However, the 
function B{f) decays quickly with / as / —?► oo. This 
guarantees that the amount of stress affecting the gene 
expression levels x does not change very rapidly. On the 
other hand, the behavior of B{f) is smooth at very small 
frequencies, i.e., B{f) —>■ i? as / —>■ 0, so that the ampli¬ 
tude of rare fluctuations of 5F remains limited. As we 
would like to understand the dynamics of x at time scales 
comparable to tMRDT, we consider the case B{f) ~ B in 
what follows. In this, the stochastic process 5F{t) in the 
right-hand side of Eq. (IA2I) has correlation properties of 
white noise: 

F{6F{t)SF{t')) = BS{t - t'), (A3) 

where 6{t — t') is the Dirac delta-function. 

3. The vicinity of a stationary point (homeostasis) 

When the fluctuations of the genotoxic forces are rela¬ 
tively weak, |5F| <C |Fb|, the fluctuations of the expres- 
some state vector are also small: x = xq + 5x^ Jcc <C xq. 
Here xq is the stationary point given by the solution of 
the equation 

ff(xo) = Fq. 

In the vicinity of this point Eq. (IA2I) reduces to 

D5x + Kdx + r^^^^xJx -!-... = SF, (A4) 

where the matrices D and K describe the relaxation ef¬ 
fects and the “interaction” between genes in the GRN. 
The matrix T^^^ encodes the leading non-linear terms in 
the expansion of the vector function g(x) in powers of 
small Sx, with all other terms corresponding to the higher 
order non-linearities being omitted from the expansion. 
Generally, such non-linear terms are not small (in partic¬ 
ular, they are not small at early times, close to the initial 


state of the expressome). However, a common feature of 
the dissipative network described by Eq. (IA4I) is the de¬ 
cay of (most) components of Sx with time, which makes 
non-linearities in Eq. (IA4I) less relevant at later stages of 
the time dynamics. Whether this feature holds for the 
GRN under consideration is a non-trivial question, which 
we shall now address. 


4. Saddle-node bifurcation and critical 
slowing-down 

Domination of a single principle component in the ex¬ 
pressome vector x{t) is not a generic prediction of the 
model (IA4I) . We would like to argue that such a domi¬ 
nance implies that GRNs under consideration are oper¬ 
ating at the critical point. This critical point is a bifurca¬ 
tion, separating regimes of stable and unstable run-away 
behavior of the expressome 5x. 

Transition from stability to instability in networks with 
network graphs not possessing specific symmetries is typ¬ 
ically associated with existence of co-dimension 1 bifur¬ 
cations [ 2 ^, [ 2 ^. Such transition is characterized by the 
loss of stability along a single orthogonal component of 
the network state vector. Gontributions from all other 
orthogonal components remain stable. This situation is 
known in the literature as saddle-node bifurcation [ 2 ^ 
and is realized when the lowest (real) eigenvalue of the 
matrix K reaches zero and then becomes negative. 

Let us denote the smallest (vanishing, but positive) 
eigenvalue of the matrix K as e and its corresponding 
left and right eigenvectors as a and b. In order to under¬ 
stand dynamics of the expressome Sx near the GRN bi¬ 
furcation, it is convenient to analyze the autocorrelation 
function F{Sx{t)Sx{t')). As shown in [sj, [s^, the lat¬ 
ter is given by the following expression (in the frequency 
domain) 

F{5x{t)5x'^ [t')) = 

= fdf- e-2"/(‘-‘')A-i(/)A(<5F(-/)<5F^(/))A*'-i(/), 

(A5) 

where A{f) = i2TrDf + K. Its behavior is determined 
by the poles of the integrand in the complex plane of 
frequency /. In turn, the latter are given by the solution 
of the equation 

det(A(/)) = 0, 

i.e., by the solution of the eigenvalue problem for the 
matrix zLi(2Tr)~^D~^K. In particular, at very large time 
separations \t — t'\ one finds that 

\ 

oA -Dh)}' 
(A6) 


F{Sx{t)Sx^ {t')) 


(a^ • Ba)b ■ b'^ 
■ Dbf 


■ exp - 
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Several important observations can be made at this point: 

1 ) as the parameter e approaches 0, amplitude of the au¬ 
tocorrelation function (jA6p grows as e“^, i.e. fluctuations 
of the expressome vector are strongly amplified both with 
age and among the population (the statistical ensemble); 

2) the autocorrelation time scale r = (a^ • Db)/e also 

grows strongly at e —>■ 0, implying that stochastic dynam¬ 
ics of the expressome becomes very slow near the point 
of bifurcation, a phenomenon known as critical slowing- 
down (see for example 3) fluctuations of the ex¬ 

pressome vector 6x mostly develop along the direction b 
of the right eigenvector of K corresponding to the van¬ 
ishing eigenvalue e, which in its turn guarantees that a 
single principal component of the expressome vector Sx 
dominates; 4) the corresponding left eigenvector a deter¬ 
mines directions of high sensitivity of the expressome to 
genotoxic stress factors SF. Other remaining eigenvalues 
of —i{27T)~^D~^K have strictly positive real parts and 
thus correspond to contributions to (IA6I) . rapidly decay¬ 
ing with time. This guarantees that only one orthogonal 
component of Sx dominates, in accordance with experi¬ 
mental observations. 

It should be noted that, when higher order derivatives 
in Eq. (lAll) are taken into account, this argument, al¬ 
beit becoming more involved, remains valid: one has to 
construct the perturbation theory in powers of small pa¬ 
rameter e, which still determines the smallest eigenvalue 
of the matrix A{f) = K + i2TTfD -|- {2nf)‘^M -|- .... Sim¬ 
ilarly to the case considered here, all higher eigenvalues 
of A{f) possess positive real parts. 

When e crosses 0 and becomes negative, the autocor¬ 
relation function (IA6I) develops unstable behavior, again 
directed along the corresponding right eigenvector b of 
K. In particular, one finds for the dispersion of the ex¬ 
pressome at large t 

E{6x{t)Sx'^{t)) « E{z^{0))b-b'^ exp (^aA^^Db) ) ’ 

where £'( 2 :^( 0 )) is the dispersion of the projection 2 = 
{Sx'^ ■ b) of the expressome on the vector 6, taken at 
the initial moment of time. The displacements of the 
expressome Sx along the other orthogonal components 
correspond to the dynamics of the gene network along 
the eigenvectors, characterizing the modes with finite and 
positive eigenvalues, and hence remain stable in time. 

Thus, the quantity of interest in the critical regime 
e —?► 0 is the projection of the expressome Sx on the right 
eigenvector b, z = {Sx'^ ■ b), which satisfies the Langevin 
equation 


^=v{z) + SE', (A8) 

where u (2) is the “velocity”, characterizing the motion 
along the 2 coordinate. Next to the origin the velocity 
v{z) Ri az, where a = — e/ {oF ■ Db). Here SF' = F/(a'^ ■ 
Db), is the stochastic force, such that E{SF'{t)SF'{t')) = 
AS{t — t'), and A = {a^ ■ Ba)/{a^ ■ Db)^ is a property of 


the fluctuations. Since the vector b is defined as eigen- 
evector and subsequently is defined up to a factor, we can 
always select its direction in such a way as to make the 
collective coordinate 2 increase with the age. It means 
that we can always limit our attention to the region 2 > 0 
only. 


Appendix B: Estimating mortality rates from the 
stochastic dynamics of an expressome 

So far we have only considered stochastic dynamics 
of the expressome state vector Sx of an organism. To 
proceed with population properties, such as mortality, 
we need to construct a statistical description for popula¬ 
tions, ensembles of organisms. Let us consider a large 
population of animals with the number of individuals 
N{t = 0) = Nq ^ 1, born at the same moment of time, 
t = 0. The development and aging of the animals in the 
population depends on the dynamics of the expression 
levels X and can be naturally described in terms of the 
probability density dN = NqP {x,t) dx, where P{x,t) is 
the fraction of individuals within the population, char¬ 
acterized by the gene expression levels from the interval 
{x, x+dx)^ and estimated at the time t. Another quantity 
of interest is the survival probability 

S(t) = j P{x,t)dx, 

defined as the fraction of the animals, N{t)/N q, surviving 
by the age t, relative to the initial size of the population 
Nq. The mortality rate can be found using the standard 
definition M{t) = —{dS/dt)/S. In this Section we shall 
explain how to estimate the probability density P(x,t)i 
the survival probability S{t), and the mortality rate M (t) 
from the stochastic dynamics of the expressome state vec¬ 
tor Sx. 


1. Estimating the probability density P{z,t) 
a. Fokker-Planck equation 


Following the arguments above we describe the process 
of aging as the slow changes in gene expression levels 
Sx over time governed by the Langevin equation (IA8I) . 
The probability P{z,t) is the solution of the associated 
Fokker-Planck equation 


where 


dP{z,t) dJjzA) 

dt dz 


(Bl) 


T / N / \ X 1 A dP{z, t) 

J {z,t) = v{z)P{z,t) - ——— 

is the probability flux along the direction 2. We seek 
the distribution function P{z, t) as the solution, corre¬ 
sponding to the initial condition P{z,t = 0) = Po{z), 
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and boundary conditions at z = 0 and large z. To ex¬ 
clude the 2 < 0 region we place a reflecting wall at z = 0, 
which corresponds to 


( dP{z,t) 
V dz 


= 0 . 

z=0 


(B2) 


The suggested form of the boundary condition means 
that the z > 0 and z < 0 regions are totally isolated 
from each other. 

We explore a Gaussian form of the initial condition in 
the form 


Po{z) = G{z) + G{-z), (B3) 


G(z) 



( 1 ^ 


characterized by the initial displacement Zq, and the dis¬ 
tribution width, (Tq. The probability distribution (IB3I) is 
an even functions and hence satisfies the boundary condi¬ 
tion for (jB2() automatically. Therefore, the normalization 
condition takes the form 


dzPo{z) = 1. 


(B4) 


In fact, in many practical situations Eq. m can be 
solved without the the boundary condition (IB2I) using a 
simpler form of the initial distribution Po{z) = G (z). In 
this case the fraction of the animals 


P_ = 


1 

2 


1-erf 







will be, of course, lost to the unphysical region z < 0. 
Therefore, the results of such a simplified calculation 
could only be trusted once the “leak” is sufficiently small, 
P_ <C 1, or, equivalently, if zq ^ cri- We distinguish be¬ 
tween the two quantitatively different types of the initial 
conditions: the “deterministic”. 


ZO » CTl, 


(B5) 


and the “diffusive”, Zq -C cti ,forms respectively. 

For sufficiently small values of z the “velocity” term in 
the Langevin equation (IA8I) can be approximated by its 
linear expansion 


V (z) = az. (B6) 

For larger z the nonlinearities set in and we employ a 
more general form 


V (z) = az(f) (z) , (B7) 

where </> (0) = 1 is chosen to match the asymptotic ex¬ 
pression (IB6I) for small z. More specifically we expect 


that the function </> (z) increases at z —> c» fast enough, 
so that that ln(z)/^(z) —?► 0 (see the discussion below). 

Biological states characterized by extremely large val¬ 
ues of z correspond to highly distorted expression pro¬ 
files and therefore we assume, that the state z = oo 
is certainly associated with the death of the individual. 
Therefore, a physically relevant choice of the boundary 
condition at large z is 


P{z —>■ oo, t) = 0. (B8) 

The Fokker-Planck equation (IBII) with the boundary con¬ 
dition in the form of Eq. (IB8I) cannot be solved exactly. 
Fortunately, this is not really needed: a large amount of 
useful information about the stochastic dynamics of the 
expressome in relation to aging can be extracted from the 
asymptotic behavior of the probability density P{z,t) in 
a few physically interesting regimes. The following anal¬ 
ysis of the solutions of Eq. m is not, of course, new 
and we loosely follow the results presented in in 

the remaining of manuscript. 


b. A toy model: absorbing wall at a large z = Z 


Let us first collect a few well known results about the 
solutions of Eq. (IBlI) for small z, where the drift velocity 
v{z) can still be used in its linear form (jB6ll . If, for 
simplicity of the argument, the initial distribution is very 
narrow, uo 0 and zq = 0 (G(z) Ri S{z)), then the 
solution satisfying the boundary condition (IB8)) is 


P{z,t) = 2G(z,0,t) = 2 


ttX (t) A 


exp 


where 


X{t)A\ ’ 
(B9) 


1 ^ 

a (z - z'e“*)^ 

V it) A 

Xit)A 


is the Green function of Eq. (IBII) . and X (t) = 
exp {2at) — 1. Therefore, at sufficiently early times, 
t < a~^, the effect of the drift term is negligible and 
the behavior of z is described by the diffusion equation 
with v{z) ~ 0. This means that at least in the beginning 
the solution is dominated by the effect of the random 
force fit) over the “deterministic” force —az 113, US- 
Over time the deterministic force is taking over. It is 
instructive to observe the relevant regimes of the solution 
using a toy model first. To emulate the effect of the non¬ 
linearity let us choose the form of the boundary condition 
at large z in the form of a totally absorbing wall, 

P(z = ±Z,t)=0, (BIO) 

suggesting that the an organism dies whenever its z 
reaches a certain threshold value Z. The solution (IB9I) 
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remains valid for sufficiently small ages, 


c. Weak non-linearity limit, general case 


t ^ In (Bll) 

before the appreciable fraction of the animals reaches the 
z ^ Z, and hence t = tis is of order of the lifespan in 
the population. Often the life expectancy tu exceeds the 
GRN instability time scale, to, ~ l/a, or 


aZ^ 


> 1 . 


(B12) 


It is shown below, that is the basic large parameter 

of our model. Whenever the ratio of the time scales is suf¬ 
ficiently large, the effects of the non-linearities are “weak” 
and a complete analytical treatment of the problem be¬ 
comes possible and yields a series of universal results. 

Technically, to meet the boundary condition (IBIOI) re¬ 
quirement, we add a loss term at the position of the ab¬ 
sorbing wall,—ry (t) [5 {z — Z) 5 {z -\- Z)], where ry (t) is 
the Lagrange multiplier, an auxiliary function introduced 
to enforce the constraint (jBlOII . The analysis of the mod¬ 
ified equation is straightforward and yields 


Let us turn our attention back to the general form of 
the non-linear drift velocity (IB7ll .At sufficiently small z 
function v[z) can be approximated by its linear expansion 
v{z) ~ az and therefore Eq. (EH) can be studied in its 
simplified form 


dp{z,t) 

dt dz 


. 1 . dP(z,t) 


0. (B17) 


We introduce a characteristic length scale z ^ Z, where 
the linear approximation to v{z) fails. The age when 
z ^ Z will correspond roughly to the lifespan of the or¬ 
ganism, since, according to our assumptions, v(z) grows 
fast enough at large 2 and the point 2; = oo is reached in 
a finite time. 

This linearized equation can be solved exactly with 
initial distribution (IB3|) and boundary condition (IB8I) : 


P{z,t) = 




■ exp < —at — 


(ze-“‘ - z^Y 

2 (t 2 


(B18) 


0 

where 


P{z,t) = <P{z,t)- 

t 

J [G{z,Z,t-t') + G{z,-Z,t-t')]r]{t')dt', (B13) 




+ 00 

{z,t) = j G{z,z',t)PQ{z')dz'. 


If the non-linearity is weak or the wall is far from cur¬ 
rent state of the system, i.e. whenever (jBI2l) holds, then 
G {Z,—Z,t — f) <C G{Z,Z,t — t'), the function rjit') 
varies very slowly with age compared to G {Z, Z,t — t'), 
and therefore 


ry (t) aZ<!> {Z, t). 

For the case of Gaussian initial distribution 


(BI4) 
we have 


^{z,t) = 


I 


•\/2^cr (t) 


a+^{t)r 

g 2,T^(t) _|_ g ^ (BI5) 


where z (t) = zoexp(at), and the effective distribution 
width increases exponentially with age 


which would be a rightful approximation for solution of 
Eq. (|Bip in time region t tig. 

In the same time, at t ^ a~^ , irrespectively of the 
specific form of the nonlinearity, one can neglect the dif¬ 
fusion and set Z\ = 0. This means that the dynamics 
of the expressome is dominated by the drift term and is, 
therefore, deterministic. Accordingly, Eq. m can be 
transformed into 




(BI9) 


One of the consequences of Ea. (IBI9p is the fact that if 
a specimen is observed with an expressome state, cor¬ 
responding to 2; 3 > a/ A/a, then the remaining life ex¬ 
pectancy, or the time-to-death, for the animal is 


00 

J viz'Y 


(B20) 


and is finite. The integral on the r.h.s. is logarithmically 
divergent for the chosen type of non-linearity EZl which 
means that 

t{z) = - In +t{z), 
a \z 


(T^ (t) =0-0-1- (f), 


(BI6) and 


where the combined quantity 



results from the initial probability density width broad¬ 
ening by diffusion at small t. 


limr (z) = const 
z—>0 


is finite. 

The solution of Ea. (jBI9p takes form P{z,t) = 
P{z,^) = G{Y)/’z{z), where ^ = t -I- t{z) and C'(^) is a 
function dependent on initial distribution and boundary 
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conditions. One way to find it is to match the solutions 
of Eqs. (IB17I) and (IB19I) at some intermediate age ti, 
lying in time region, where both equations can be used. 
Indeed, since normally a~^ tjs > we can use any t\ 
such that 


a ^ < tis (B21) 


as the matching point. The age range corresponds to 
A/a zi < Z interval. As the result, we obtain a 
closed form solution, applicable for deterministic initial 
condition (IB5I) and boundary condition (IB8I) for all ages 
t > 1/a 


P{z,t) 


Here 


ak (z) f [k (z) e 

,_ _exp < —at - 

\f^av (z) 2 (t2 


fof I 

(B22) 


Z = lime“^(^) - Z, (B23) 

0—>0 

and 

fc (z) = Zexp [—ar (z)] (B24) 

is auxiliary function with known asymptotic forms for 
the advanced (r (z) —>■ 0) 


k{z) ^ Z at z ^ oo. 

(B25) 

and the early (r (z) —>■ oo) 


fc (z) —?► 0 at z —>■ 0 

(B26) 

ages. 

Let us turn to a more specific case of an arbitrary poly¬ 
nomial non-linearity 


(B27) 


with /3 > 0. The length scale here characterizes the non¬ 
linearity strength and hence large values of Z correspond 
to the weak non-linearity limit (jB6() . In accordance with 
the earlier definitions, it can be shown that Z equals Z, 
and the remaining lifespan is 




and the probability distribution function takes the form 

of Eq. (IB22|) with with fe (z) = z/ [l -I- (z/Z)^] 

The absorbing wall at z = Z situation, considered in 
the previous section, formally corresponds to /3 = oo. 
Interestingly, the solution is still valid, k{z) = z every¬ 
where, if not immediately close to the wall, i.e. Eq. (IB22I) 
matches the exact expression (IBI3I) for all z < Z. Hence 


we can conclude that the derived form of the probabil¬ 
ity density is universal, holds for an arbitrary form of 
the non-linearity, allowing for a finite lifespan. A specific 
form of the non-linearity is not important, all the infor¬ 
mation about the non-linear terms can be compressed 
into the length scale Z, being the single important pa¬ 
rameter defining the lifespan and the form of the distri¬ 
bution function in the population. 

Let us finish the section by computing the fraction of 
the animals, surviving by the age t. According to the 
definition. 


+ CX) 

S {t) = j dzP (z, t) 


^ r I \i. (^\ 




dk exp < —at — 


[k (z) e-“* - zo]^ 


2 a 2 


y2 

: J dye~y^, 


yi 

where, for at ^ 1, 

_ Ze~°‘* - zo _ -zq 

That would lead us to expression 


-S' (0 = ^ ^ erf 


.V2cr. 


- erf 


- ZQ 

V2a 


(B28) 


We also note, that 


+ 00 

J dzP (z, t) < I, 


since the animals die and hence disappear at the infinity. 


2. Estimating the first passage time, the survival 
probabilities, and the average lifespan 

In the case under consideration, a more interesting 
quantity to calculate and analyze is not the probability 
density P{z,t) itself, but a so-called first passage time 
distribution function PFPT(t) (H, The reason is 

that the probability density P(z, t) represents a sum over 
all possible stochastic trajectories z(t), including those 
which pass the effective threshold z = Z and then return 
back. On the other hand, by definition, Pfpt(T) is a 
probability for a given stochastic trajectory z = z{t) to 
pass the threshold z = Z at t = T for the very first time. 
As such, it does not take into account trajectories which 
then re-enter the allowed domain of z. 
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When the “deterministic” form of the initial condition 
is chosen, the first passage time distribution function is 
given by 


PFPT(t) = —S {t) = 


1 a 2aZ'e^“* 

aZ2 

V ttA [X (t)]3/2 

X(t)A_ 


(B29) 


The function (IB29I) grows exponentially 
a~^, possesses a distinct maximum at t 


In 


( 

V ^ 


1 / 2 - 


with its width 6t 


at small t < 

max — ils — 

a~^ near the 


maximum and finally decays as e““* at t ^ tmax, see Fig. 
O Computing the first passage time distribution function 
also allows one to easily estimate the survival probability 
S{t), using the prescription Ppprit) = —dS{t)/dt. From 
Eq. (IB29I) we find 


S{t) = erf 


I aZ^ 

2AX(t) 


(B30) 


In the regime a ^ < t ^ tis , also known as Suzuki scaling 
[13, [13 ,equation for Ppprit) reduces to 



Figure 5: The first passage time distribution function 
PppTit) (in units of 2a) and the survival probability 
S{t) as a function of time 2at; = 1000. 


PFPT(t) 



(B31) 

If the non-linearity is small, the center of the originally 
narrow distribution in the z-space follows the well defined 
trajectory z (t) = zq exp (at). At t = to = a~^ In (Zj zq) 
it achieves the absorbing wall, which means that the aver¬ 
age lifespan in the population described by the distribu¬ 
tion is also tis Pi to- This result follows immediately from 
(IB28I) and the exact definition of the average lifespan 


OO 



0 


From Eq. (|B3ip we find 


(B32) 


S(t) ~ erf 



(B33) 


expression precise up to terms, exponentially small at 
^ 1, i.e., when a species is considered long-lived 
compared to the characteristic time scales of molecular 
dynamics in the cell, see the next Subsection. The func¬ 
tion (|B33p . initially exponentially close to 1, then falls off 
quickly (as we shall see below, this fall-off corresponds to 
an exponentially growing mortality rate) and finally ap¬ 
proaches zero as exp(—at) at t ^ tig. 


3. Behavior of the mortality rate 


As usual, the mortality rate is defined as M{t) = 
—S~^ ■ dS/dt, so that one has precisely 


Mit) = (B34) 

Once again, for the very narrow initial distribution (IB8p . 
CTo —>■ 0 and Zq = Q (G{z) ~ d{z)), we find, that 


M(t) 


2aZe2“* 

V 7rZ\ [X 


■ aZ^ 

~ X (t) A _ 


X 


X 



' aZ^ 

2X{t)A 


-1 


(B35) 


For a more realistic determenistic form of the probability 
distribution at the origin lB5l the solution is still possible, 
and for sufficiently advanced ages, t ^ a“^, we obtain 


M{t) 


Za 

<j 



(Ze-“‘ - zo)' 

2 a 2 


— at 


X 


X 


11 -I- erf 


Ze-“* - zo 



(B36) 


once the non-linearity is sufficiently weak, as commanded 
by the condition (IB12I) . 
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a. Vanishing mortality at very early ages t a 

At sufficiently early ages t <C q;“^ behavior of z is 
dominated by the effect of the random force f{t) over 
the “deterministic” force v{z) ^t\. Accordingly, the 
state variable and its fluctuations grow exponentially 
with age. This conclusion is confirmed by the analysis 
of senescence-associated transcriptional dynamics. It is 
worth emphasizing that the characteristic time scale of 
the stochastic growth of the expressome is ~ q;“^ . It thus 
depends only on the morphological properties of the GRN 
and does not depend on the genotoxic stress amplitude 
A at all. 

At very early ages t a~^ the first passage time dis¬ 
tribution function is exponentially strongly suppressed: 

In this regime the survival probability S(t) remains ex¬ 
ponentially close to 1: 


first passage time probability distribution Pfpt reaches 
a maximum at tmax = In A = tis, where A = 
and its half-width near the the maximum is of the order 
a~^, as in Fig. O By the time t ~ tmax the argument 
of the error function in Eq. (|B33I) is already sufficiently 
small for it to be expanded. One thus has 

S{t) ~ -^Ae““‘-^A^exp(—3at) -I-_ 

VTT 3V7r 

On the other hand, one finds for PFPT(t) in the same 
regime 

PfptW = 
so that approximately 
M(<) « 

Introducing t = tmax + St one finally finds at small St < 
(2a)-i 


S{t) 





Therefore, the mortality rate remains exponentially close 
to 0, and one approximately has M{t) r; PFPT(t)- In our 
opinion, this behavior might explain relatively small mor¬ 
tality, often observed in cohorts at very early ages. We 
note that the behavior of the mortality rate in this regime 
is not universal as it strongly depends on the amplitude 
A of stochastic genotoxic stress. 


M{t) « _ o.61ae“^*. (B38) 

Thus, the Gompertz law with the universal exponent a~^ 
is realized at times t^ax—(2 q;)“^ ^ t ^ imax+(2Q;)“^. We 
again emphasize that the Gompertz aging rate a depends 
only on the properties of the GRN under consideration 
and does not depend on stress. In fact, the very same 
conclusion holds for the Gaussian initial distribution case 

(IB361) . 


b. Gompertz behavior at later ages t ^ a 


As was explained above, as the age increases, the effect 
of the stochastic force SF in Eq. (IA4I) quickly becomes 
negligible in comparison to the effect of the “determin¬ 
istic” force —KSx; the stochastic regime is changed by 
the drift motion at t ~ a~^. Correspondingly, the be¬ 
havior (IB37|) quickly changes when t > ct~^, yielding a 
Gompertz-like dependence of the mortality rate on age: 
M [t) ft! Mo exp (qt), where the exponent 


lit) = w 


aZ^X (f) 
A 2 (t) A 



> a 


is found from Eq. (IB35|) and is much greater than the 
GNR instability parameter a. 


d. Mortality deceleration at late ages 

Finally, at t ^ t^ax both the first passage time distri¬ 
bution function (IB29I) and the survival probability (IB30I) 
decay exponentially. This corresponds to an asymptot¬ 
ically constant behavior of the mortality rate in both 
cases, (jB35ll and (jB36ll : 

M{t) ft a, (B39) 

see Fig. [HI For late t, corrections to the leading constant 
term become exponentially small with time. The exis¬ 
tence of mortality plateau in models of the type (|A4p is 
known in literature . 


4. Concluding remarks on aging regimes 
phenomenology 


c. Universal Gompertz behavior in the regime of Suzuki 
scaling 

At even later times of the order t ^ tig, when the 
regime of Suzuki scaling is realized, we also find a uni¬ 
versal Gompertz behavior. As was explained above, the 


a. The unstable case, Gompertz limit, aZ'^ jA 1 

This case includes situations with a very weak non¬ 
linearity in the Fokker-Planck equation (IBlI) and corre¬ 
sponds to a very high toxicity threshold Z as compared to 
the characteristic amplitude of fluctuations of the stress 
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Figure 6: The mortality rate M{t) for the case of the 
“diffusive” initial condition; aZ'^ jA = 1000. The 
maximum of PppTit) is reached at 2at Ri 6.9. 



Figure 7: The mortality rate as a function of age from 
the very large cohorts study, medflies, males and 
females from [2l|. 


factors A. As was explained above, we expect a species 
with such a GRN to be relatively long-lived, as the aver¬ 
age life expectancy (IBllI) is logarithmically larger than 
the inverse Gompertz exponent a~^. Four very distinct 
regimes of mortality exist for such species: (a) the regime 
of small mortality at early ages, (b) a Gompertz regime 


with the Gompertz exponent ^ o?Z'^ jA ^ a, (c) a Gom¬ 
pertz regime with a universal Gompertz exponent a. and 
finally (d) a regime of constant mortality M^o = a. On 
the other hand, if indeed Z ^ ctoj the initial Gompertz 
exponent in regime (b) is very large, and a population of 
moderate size simply dies out before the regime of mor¬ 
tality slowing-down is observed. 

h. The unstable case, short lifespan limit, aZ'^ jA > 1 

There are again four distinct regimes (a-d) of aging 
discussed above, with that difference that the regime 
(b) is hard to observe or unobservable, so a Gompertz 
regime (c) with the exponent ~ a is observed. Glearly, 
the regimes (c) and (d) are closely connected to each 
other in the sense, described above, and there is a strong 
correlation between the observed value of the Gompertz 
exponent and the mortality rate Moo at late ages. This, 
as it seems, is what we observe e.g. in experiments with 
extremely large cohorts study of medflies, males and fe¬ 
males from [^, see Fig. [T] A relatively short aver¬ 
age lifespan of such species is comparable to the value of 
Gompertz exponent. 

c. The stable case, a < 0 

In this case, arguably the most interesting to consider, 
the effective potential V{z) of the Langevin equation pos¬ 
sesses a local minimum, separated from the roll-down to 
2 —>■ cxD by a potential barrier, see Fig. [TJ3 of the main 
text. There exists a single regime of aging realized for 
the GRN of this type. It is characterized by a constant 
mortality rate Moo, its value being determined by the ex¬ 
ponentially small Kramers probability of passage through 
the barrier If a GRN with such properties indeed ex¬ 
ists in nature, it would correspond to a species with neg¬ 
ligible senescence. We leave the in-depth study of this 
regime for a future work. 

Appendix C: Gene annotation 

To analyze the vector 6, obtained from the transcrip- 
tome of D. melanogaster [l9|, we have performed a gene 
set enrichment analysis for the components of vector b us¬ 
ing “Parent-Child-Intersection” GO enrichment analysis 
procedure described in (^ . This method considers genes 
in the background if they are present in all parent terms 
of the term of interest. Since the values in vector b can 
be both positive and negative, lists of positive and neg¬ 
ative leading components were analyzed separately. We 
studied the dependence of significantly enriched GO cate¬ 
gories on the number of the leading components selected. 
The list of enriched GO categories appeared to be stable 
in a broad range of number of the leading components, so 
200 positive and 200 negative leading components from 
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GO-stability regions were analyzed.We employed an ad¬ 
justed p-value (Benjamini-Hochberg) cutoff of 0.05. 

We have performed KEGG pathway enrichment anal¬ 
ysis of the human orthologs of the genes corresponding 
to the most significant components of vector b, computed 


from the transcriptome of D. melanogaster [l^ . A one¬ 
sided Fischer exact test was performed, after which an 
adjusted p-value (Benjamini-Hochb erg ) cutoff of 0.05 was 
employed. We used OrthoDB v.7 to map fly genes 
to corresponding human orthologs. 
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